//==============================================================================
// SeqSampler.cc
// ----------------------------------------------------------------------------
// 
//==============================================================================
//.......1.........2.........3.........4.........5.........6.........7.........8
//

#include "SeqSampler.hh"

RandomNumberGenerator RNG;


bool
SeqSampler::OnlyBases(Sequence S, int idx,int fsize){
  for(int b = 0;b < fsize;b++){
	 char letter = toupper(S[idx+b]);
	 if ((letter == 'N') || (letter == 'X')) return(false);
  }
  return(true);
}

//------------------------------------------------------
// * FastFragmentSeq
//------------------------------------------------------
// Attempt to make a faster FragmentSeq with restricted options.  
//
//
void 
SeqSampler::FastFragmentSeq(int FLength,Sequence &SeqIn,SequenceVector &SVOut){

  int seqSize = SeqIn.size();
  int numFrags = seqSize-FLength+2;

  // Create space for samples
  SVOut.resize(numFrags);
  // For each sample...
  int fStart = 0;
  for(int i = 0;i< numFrags;i++){

	 // Skip samples that contain N or X. 
	 if (OnlyBases(SeqIn,fStart,FLength)){
		SVOut[i].resize(FLength);
		// Copy the sample from the input sequence
		for(int b = 0;b< FLength;b++){
		  SVOut[i][b]=SeqIn[fStart+b];
		}
	 }
	 fStart++; // Move start over one base;

	 //SVOut[i].mDescription = SeqIn.mDescription+itoa(i); // KJD leave description out for now. 
  }
}

//------------------------------------------------------
// * FragmentSeq
//------------------------------------------------------
// Completely breaks a sequence up into fragments, which 
// overlap by FOverlap.   FOverlap may be zero, but may not be negative. 
// 
void
SeqSampler::FragmentSeq(int FLength,int FOverlap,Sequence &SeqIn,
								SequenceVector &SVOut)
{
  if (FOverlap > FLength) return;
  Sequence tempSeq;
  long fStart=0,fEnd=FLength-1;
  int SeqSize = SeqIn.size();
  do{
	 tempSeq.clear();
	 SeqIn.GetSubSequence(tempSeq,fStart,fEnd);
	 //SVOut.resize(SVOut.size()+1);
	 //SVOut[SVOut.size()-1]=tempSeq;
	 SVOut.push_back(tempSeq);
	 fStart = fEnd-FOverlap+1;
	 fEnd = fStart+FLength-1;
  }while(fEnd<SeqSize);
}

//------------------------------------------------------
// * FragmentSeq
//------------------------------------------------------
//
void
SeqSampler::FragmentSeq(int FLength,int FOverlap,Sequence &SeqIn,
								SequenceVector &SVOut,
								vector<long> &StartLocs,vector<long>&EndLocs)
{
  if (FOverlap > FLength) return;
  Sequence tempSeq;
  long fStart=0,fEnd=FLength-1;
  int SeqSize = SeqIn.size();
  do{
	 tempSeq.clear();
	 StartLocs.push_back(fStart);
	 EndLocs.push_back(fEnd);
	 SeqIn.GetSubSequence(tempSeq,fStart,fEnd);
	 SVOut.push_back(tempSeq);
	 fStart = fEnd-FOverlap+1;
	 fEnd = fStart+FLength-1;
  }while(fEnd<SeqSize);
}


//------------------------------------------------------
// * RndFragmentSeq
//------------------------------------------------------
//
void
SeqSampler::RndFragmentSeq(int FLength,int FOverlap,Sequence &SeqIn,
									SequenceVector &SVOut,
									vector<long> &StartLocs,vector<long>&EndLocs,
									vector<long> &Overlaps)
{
  if (FOverlap > FLength) return;
  Sequence tempSeq;
  long fStart=0,fEnd=FLength;
  int SeqSize = SeqIn.size();
  do{
	 cout<<"fStart="<<fStart<<" fEnd="<<fEnd<<" seqSize="<<SeqSize<<flush;
	 tempSeq.clear();
	 StartLocs.push_back(fStart);
	 EndLocs.push_back(fEnd);
	 SeqIn.GetSubSequence(tempSeq,fStart,fEnd);
	 SVOut.push_back(tempSeq);
	 long tempOverlap = (long)RNG.Gaussian(FOverlap,0.3*FOverlap);
	 if (tempOverlap < 100) tempOverlap = 100;
	 if (tempOverlap > (0.5*FLength)) tempOverlap = (long)(0.5*FLength);
	 Overlaps.push_back(tempOverlap);
	 cout<<"overlap = "<<tempOverlap<<endl;
	 fStart = fEnd-tempOverlap-1;
	 fEnd = fStart+FLength-1;
  }while(fEnd<SeqSize);
}


//---------------------------------------------------
// * SampleSeq
//---------------------------------------------------
// Generates potentially non-overlapping samples of a 
// sequence. 
//
void
SeqSampler::SampleSeq(int NumSamples,int SampleSize,
							 Sequence &SeqIn,SequenceVector &SamplesOut)
{
  long SampleSpacing = SeqIn.size()/NumSamples;
  for(long sNum = 0;sNum < NumSamples;sNum++){
	 long start = (long)(sNum*SampleSpacing);
	 long end = start+SampleSize - 1;
	 SamplesOut.AddSubSequence(SeqIn,start,end);
  }
}



//------------------------------------------------------
// * AddReverseComplements
//------------------------------------------------------
//
void
SeqSampler::AddReverseComplements(SequenceVector &SeqVec)
{
  Sequence temp;
  int OriginalSize = SeqVec.size();
  for(int i = 0;i< OriginalSize;i++){
	 temp = SeqVec[i];
	 temp.ReverseComplement();
	 temp.mName += "(-)";
	 SeqVec.push_back(temp);
  }
}

//------------------------------------------------------
// * CreateReverse
//------------------------------------------------------
//
void
SeqSampler::CreateReverse(SequenceVector &SeqsIn,
								  SequenceVector &SeqsOut)
{
  Sequence temp;
  SeqsOut.clear();
  SeqsOut.resize(SeqsIn.size());
  for(int i = 0;i< SeqsIn.size();i++){
	 temp = SeqsIn[i];
	 temp.Reverse();
	 temp.mName +="(r)";
	 SeqsIn.push_back(temp);
  }
}

//------------------------------------------------------
// * CreateReverseComplement
//------------------------------------------------------
//
void
SeqSampler::CreateReverseComplement(SequenceVector &SeqsIn,
  SequenceVector &SeqsOut)
{
  Sequence temp;
  SeqsOut.clear();
  SeqsOut.resize(SeqsIn.size());
  for(int i = 0;i< SeqsIn.size();i++){
	 temp = SeqsIn[i];
	 temp.ReverseComplement();
	 temp.mName +="(-)";
	 SeqsIn.push_back(temp);
  }
}


//------------------------------------------------------
// * CreateComplement
//------------------------------------------------------
//
void
SeqSampler::CreateComplement(SequenceVector &SeqsIn,
												SequenceVector &SeqsOut)
{
  Sequence temp;
  SeqsOut.clear();
  SeqsOut.resize(SeqsIn.size());
  for(int i = 0;i< SeqsIn.size();i++){
	 temp = SeqsIn[i];
	 temp.Complement();
	 temp.mName +="(c)";
	 SeqsIn.push_back(temp);
  }
}


//------------------------------------------------------
// * AddAllSenses
//------------------------------------------------------
//
void
SeqSampler::AddAllSenses(SequenceVector &SeqVec)
{
  Sequence temp;
  int OriginalSize = SeqVec.size();
  for(int i = 0;i< OriginalSize;i++){
	 temp = SeqVec[i];
	 temp.ReverseComplement();
	 temp.mName += "(-)";
	 SeqVec.push_back(temp);

	 temp = SeqVec[i];
	 temp.Reverse();
	 temp.mName += "(r)";
	 SeqVec.push_back(temp);

	 temp = SeqVec[i];
	 temp.Complement();
	 temp.mName += "(c)";
	 SeqVec.push_back(temp);	 
  }
}


//---------------------------------------------------
// * GetRandomSample
//---------------------------------------------------
// Grabs a random sample of the given size from a sequence 
//
void
SeqSampler::GetRandomSample(int sSize,Sequence &Source,Sequence &SampleOut)
{
  RandomNumberGenerator RNG;
  long maxStart = (long)(Source.size()-(long)sSize);

  // Start of a read is uniformly distributed over the target molecule
  if (maxStart < 0){
	 cout<<"maxStart < 0!!!!!"<<maxStart<<endl;
	 cout<<"Source.size()="<<Source.size()<<endl;
	 cout<<"sSize="<<sSize<<endl;
	 exit(1);
  }

  long sstart = (int)RNG.Uniform(0,maxStart);
  long send = sstart+sSize-1;

  if (sstart > maxStart) sstart = maxStart;
  if (sstart < 0) sstart = 0;
  if (send > (Source.size()-1)) send = Source.size() -1;

  Source.GetSubSequence(SampleOut,sstart,send);

  if (SampleOut.size()==0){
	 cout<<"1start="<<sstart<<endl;
	 cout<<"1end="<<send<<endl;
	 cout<<"maxStart="<<maxStart<<endl;
	 cout<<"Source.size="<<Source.size()<<endl;
  }

}







